ohi logo
OHI-Northeast | OHI Science | Citation policy


#Data Source

Downloaded: December 14, 2018 (emailed by Jeffrey Vieser at NMFS)

Description: Records of Bmsy and Fmsy estimates from stock assessments conducted in the greater Northeast region

Time range: 2004 - 2018

Format: Tabular


1 Data wrangling

1.4 Span plot

Look at stock assessment data by stock over time to see what we are working with.

What species have multiple stocks? These ones will give us a tough time when matching up catch to the stocks.

## # A tibble: 20 x 3
##    stock                             species        location               
##    <chr>                             <chr>          <chr>                  
##  1 Atlantic cod - Georges Bank       Atlantic cod   Georges Bank           
##  2 Atlantic cod - Gulf of Maine      Atlantic cod   Gulf of Maine          
##  3 Haddock - Georges Bank            Haddock        Georges Bank           
##  4 Haddock - Gulf of Maine           Haddock        Gulf of Maine          
##  5 Windowpane - Gulf of Maine / Geo… Windowpane     Gulf of Maine / George…
##  6 Windowpane - Southern New Englan… Windowpane     Southern New England /…
##  7 Winter flounder - Gulf of Maine   Winter flound… Gulf of Maine          
##  8 Winter flounder - Southern New E… Winter flound… Southern New England /…
##  9 Yellowtail flounder - Cape Cod /… Yellowtail fl… Cape Cod / Gulf of Mai…
## 10 Yellowtail flounder - Georges Ba… Yellowtail fl… Georges Bank           
## 11 Yellowtail flounder - Southern N… Yellowtail fl… Southern New England /…
## 12 Winter flounder - Georges Bank    Winter flound… Georges Bank           
## 13 Red hake - Gulf of Maine / North… Red hake       Gulf of Maine / Northe…
## 14 Red hake - Southern Georges Bank… Red hake       Southern Georges Bank …
## 15 Silver hake - Gulf of Maine / No… Silver hake    Gulf of Maine / Northe…
## 16 Silver hake - Southern Georges B… Silver hake    Southern Georges Bank …
## 17 Atlantic cod - Eastern Georges B… Atlantic cod   Eastern Georges Bank   
## 18 Haddock - Eastern Georges Bank    Haddock        Eastern Georges Bank   
## 19 Goosefish - Gulf of Maine / Nort… Goosefish      Gulf of Maine / Northe…
## 20 Goosefish - Southern Georges Ban… Goosefish      Southern Georges Bank …

Looks like there are just 8.

3 Calculate stock scores

Convert these value to stock scores.

#grab just B/Bmsy data
  nmfs_b_bmsy <- data_gf %>%
    filter(indicator == "b_bmsy") %>%
    select(stock, year, value)

#grab just F/Fmsy data  
  nmfs_f_fmsy <- data_gf %>%
    filter(indicator == "f_fmsy") %>%
    select(stock, year, value)
  ### 
  b_bmsy_underexploit_penalty <- 0.25
  b_bmsy_underexploit_thresh  <- 3.00
  f_fmsy_underfishing_penalty <- 0.25
  f_fmsy_overfishing_thresh   <- 2.00
  
  ### Apply rolling mean to F/Fmsy
  ## Why do we do this? because B is a less sensitive metric (relies of biological processes) and F can fluctuate pretty easily because it is really just a mgmt decision.
  nmfs_f_fmsy <- nmfs_f_fmsy %>%
    arrange(stock, year) %>%
    group_by(stock) %>%
    mutate(f_fmsy_rollmean = zoo::rollmean(value, k = 3, align = 'right', fill = NA)) %>%
    ungroup() %>%
    select(-value) %>%
    rename(value = f_fmsy_rollmean)
  
  stock_status_layers <- nmfs_b_bmsy %>%
    full_join(nmfs_f_fmsy) %>%
    spread(indicator, value) 
  
########################################################.
##### run each fishery through the Kobe plot calcs #####
########################################################.
### * ram_b_bmsy, ram_f_fmsy
  
  
### Function for converting B/Bmsy values into a 0 - 1 score
  rescale_bprime_crit <- function(fish_stat_df,
                                  bmax, bmax_val) {
    
    ###using NOAA's limits here
    overfished_th  <- 0.8
    ### 
    underfished_th <- 1.2
    
    bmax_adj <- (bmax - underfished_th) / (1 - bmax_val) + underfished_th
    ### this is used to create a 'virtual' B/Bmsy max where score drops
    ### to zero.  If bmax_val == 0, this is bmax; if bmax_val > 0, bmax_adj
    ### extends beyond bmax, to create a gradient where bmax_val occurs at bmax.
    
    fish_stat_df <- fish_stat_df %>%
      # group_by(stock) %>% ### grouping by stock will set b_max by max per stock, instead of max overall
      mutate(b_max     = max(b_bmsy, na.rm = TRUE)) %>%
      ungroup() %>%
      mutate(bPrime = NA,
             bPrime = ifelse(b_bmsy < overfished_th,  ### overfished stock
                             b_bmsy / overfished_th,
                             bPrime),
             bPrime = ifelse(b_bmsy >= overfished_th & b_bmsy < underfished_th,
                             1,                       ### appropriately fished stock
                             bPrime),
             bPrime = ifelse(b_bmsy >= underfished_th,
                             (bmax_adj - b_bmsy) / (bmax_adj - underfished_th), ### underfished stock
                             bPrime),
             bPrime = ifelse(bPrime < 0, 0, bPrime))
    return(fish_stat_df)
  }
  
  
  ### Function to create vertical gradient based on distance from
  ### ideal F/Fmsy value to actual F/Fmsy value
  f_gradient <- function(f, over_f, under_f, fmax, fmin_val) {
    x <- ifelse(f < over_f & f > under_f, 1, NA)
    x <- ifelse(f <= under_f, (f * (1 - fmin_val) / under_f + fmin_val), x)
    x <- ifelse(f >= over_f,  (fmax - f) / (fmax - over_f), x)
    x <- ifelse(f > fmax, NA, x)
    return(x)
  }
  
  ### Function to convert F/Fmsy values into 0 - 1 score
  rescale_fprime_crit <- function(fish_stat_df,
                                  fmax, fmin_val) {
    
    ### params - taken from BC but changed Bcrit to 0.5 instead of 0.4:
    Bcrit <- 0.5; overfished_th <- 0.8
    ### underfishing_th is set to the idea "1/3 for the birds":
    underfishing_th <- 0.66; overfishing_th  <- 1.2
    
    bcritslope = 1 / (overfished_th - Bcrit)
    ### connecting from (Bcrit, 0) to (overfished_th, 1)
    
    fish_stat_df <- fish_stat_df %>%
      mutate(fPrime = ifelse(b_bmsy < overfished_th & f_fmsy < fmax,
                             f_gradient(f_fmsy + (overfished_th - b_bmsy) * bcritslope,
                                        over_f = overfishing_th,
                                        under_f = underfishing_th,
                                        fmax = fmax,
                                        fmin_val = fmin_val),
                             NA),
             fPrime = ifelse(b_bmsy >= overfished_th & f_fmsy < fmax,
                             f_gradient(f_fmsy,
                                        over_f = overfishing_th,
                                        under_f = underfishing_th,
                                        fmax = fmax,
                                        fmin_val = fmin_val),
                             fPrime),
             fPrime = ifelse(is.na(fPrime), 0, fPrime), ### fill zeros everywhere unscored
             fPrime = ifelse(is.na(f_fmsy), NA, fPrime) ### but if no f_fmsy, reset to NA
      )
    return(fish_stat_df)
  }
  
  stock_status_df <- stock_status_layers %>%
    rescale_bprime_crit(bmax     = b_bmsy_underexploit_thresh,
                        bmax_val = b_bmsy_underexploit_penalty) %>%
    rescale_fprime_crit(fmax     = f_fmsy_overfishing_thresh,
                        fmin_val = f_fmsy_underfishing_penalty) %>%
    mutate(x_prod = ifelse(!is.na(fPrime), (fPrime * bPrime), bPrime),
           basis  = case_when(
             !is.na(fPrime) & !is.na(bPrime) ~ 'F_Fmsy, B_Bmsy',
             is.na(fPrime)  & !is.na(bPrime) ~ 'B_Bmsy only',
             is.na(bPrime)  & !is.na(fPrime) ~ 'F_Fmsy only'
           )) %>%
    dplyr::select(year, stock,
                  score = x_prod,
                  basis,
                  bPrime, fPrime,
                  b_bmsy, f_fmsy) 

Take the average scores for the 5 species with sub-stocks. Just doing a group_by with species and averaging scores will do this for all species (but only change the values for these sub-species).

Remove stocks with just F_Fmsy - these stocks have no stock scores since we don’t have a way to get a score from just F/Fmsy. Also roll all values forward to 2017.

4 Visualize

kobe_plot <- function(sp){
  
 ss_df <- stock_scores %>%
    filter(stock == sp) %>%
    arrange(year) %>%
    mutate(last_bbmsy = last(b_bmsy),
           last_ffmsy = last(f_fmsy),
           last_datayear = last(year)) %>%
   ungroup()

generate_kobe_df <- function(f_fmsy_max = 2.5,
                             b_bmsy_max = 3.0,
                             reso       = 0.01,
                             bmax_val   = 0,
                             fmin_val   = 0,
                             weighting_b = 1) {

  kobe_raw <- data.frame(stock  = 1,
                     f_fmsy = rep(seq(0, f_fmsy_max, reso), each  = round(b_bmsy_max/reso) + 1),
                     b_bmsy = rep(seq(0, b_bmsy_max, reso), times = round(f_fmsy_max/reso) + 1))

  kobe <- kobe_raw %>%
    rescale_bprime_crit(bmax = b_bmsy_underexploit_thresh,
                        bmax_val = bmax_val) %>%
    rescale_fprime_crit(fmax = f_fmsy_overfishing_thresh,
                        fmin_val = fmin_val) %>%
    mutate(x_geom  = (fPrime * bPrime),
           x_arith = (fPrime + bPrime) / 2)

  return(kobe)
}

bbmsy_lim <- max(round(max(ss_df$b_bmsy, na.rm = TRUE) + .1, 1), 3)
ffmsy_lim <- max(round(max(ss_df$f_fmsy, na.rm = TRUE) + .1, 1), 2.5)
  
kobe_df <- generate_kobe_df(f_fmsy_max = ffmsy_lim,
                           b_bmsy_max = bbmsy_lim,
                           bmax_val = .25,
                           fmin_val = .25)


kobe_stock_plot <- ggplot(data = kobe_df, aes(x = b_bmsy, y = f_fmsy)) +
    theme_bw() +
    geom_raster(alpha = .8, aes(fill = x_geom)) +
    scale_fill_distiller(palette = 'RdYlGn', direction = 1) +
    labs(title = as.character(sp),
         x = 'B/Bmsy',
         y = 'F/Fmsy',
         fill = "Stock score") +
    geom_path(data = ss_df, 
              show.legend = FALSE,
              aes(x = b_bmsy, y = f_fmsy, group = sp),
              color = 'grey30') +
    geom_point(data = ss_df, 
               show.legend = FALSE,
              aes(x = last_bbmsy, y = last_ffmsy)) +
    geom_text(data = ss_df %>%
                mutate(year = ifelse(year/5 == round(year/5) | year == last_datayear, year, NA)), 
              aes(x = b_bmsy, y = f_fmsy, label = year), 
              hjust = 0, nudge_x = .05, size = 2)

return(kobe_stock_plot)
}

Apply function to all species in the stock_scores data frame.